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Abstract. Frequency conversion (FC) and parametric down-conversion (PDC) 
are among the most widely used nonlinear processes for the implementation of 
quantum optical experiments. Parametric down- conversion enables the efficient 
creation of quantum states ranging from photon-number states over squeezers to 
EPR-states. Frequency conversion gives rise to technologies enabling efficient 
atom-photon coupling, ultrafast pulse gates and enhanced detection schemes. 
Flowever, despite their widespread deployment, their theoretical treatment 
remains challenging. Especially the multi-photon components in the high gain 
regime, as well as the explicit time-dependence of the involved Hamiltonians 
hamper an efficient theoretical description of these nonlinear optical processes. 

In this paper we investigate these effects and put forward two models which 
enable a full description of FC and PDC in the high gain regime. We present a 
rigorous numerical model relying on the solution of coupled integro-differential 
equations which covers the complete dynamics of the process. And, as an 
alternative, we develop a simplified model that, at the expense of neglecting time- 
ordering effects, enables an analytical solution which approximates the correct 
solution with high fidelity in a broad parameter range. 

The developed fundamental understanding of frequency conversion and 
parametric down- conversion gives valuable insights into the quantum properties 
of the processes, extends the current theoretical descriptions, and simplifies 
considerably the engineering process for future quantum information applications 
using FC and PDC. 
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1. Introduction 

A fundamental building block of quantum information and quantum communication 
applications are nonlinear optical processes. In experimental implementations 
of photonic quantum systems parametric down-conversion (PDC) and frequency 
conversion (FC) are omnipresent. Parametric down-conversion enables the generation 
of various quantum states ranging from single photons [U [2j [3l HI [5] over entangled 
photon-pairs [5J [71 [5] up to squeezers [HI [TU] and EPR-states [TTJ [T2]. Frequency 
conversion is applied for frequency translations between different wavelengths [131 HH 
I15j which enables interfaces between quantum systems, in particular atom-photon 
coupling [ini E], quantum pulse gates [TSl [HI [50], and efficient quantum state 
detection [HJ [22l [H [H [25] . 

Their deployment in quantum enhanced applications requires a detailed 
theoretical understanding of the corresponding nonlinear interactions. A variety of 
models have been developed for PDC [Ml[2a[Ml[Ml[l2l[3Q]andFC[18l[19l[20l[l5l[31]. 
They vary from straightforward perturbation approaches to much more rigorous 
treatments. The crucial issue in these derivations is firstly the fact that multi-photon 
effects have to be considered during the interaction, and secondly the problem that 
the involved electric field operators and consequently Hamiltonians do not commute 
in time. In this paper we address these issues and build two theoretical models for 
FC and PDC: a rigorous numerical model based on the theoretical framework of 
Mikhail I. Kolobov [3^, and a simplified analytical approach. Both models take into 
account higher-order photon number effects and are hence suitable to describe FC and 
PDC in the high gain regime. We analyse their performance and the quality of their 
predictions over a broad parameter range which enables us to suggest when a simple 
analytic modelling of the processes is sufficient and when the rigorous approach has 
to be applied. 

The paper is structured into two main parts: In sections [2] to [9] we study frequency 
conversion. Our investigation of this process is divided into eight sub chapters: After a 
short description of the basic principles of frequency conversion in section [21 section [3l 
discusses the Hamiltonian of the process. The general properties of the conversion are 
outlined in section [31 In section [SI we introduce the no time-ordering approximation 
which forms the basis for the analytical solution of the frequency conversion process 
and elaborate on its implications. In scction[6lwe derive the analytic solution excluding 
time-ordering effects. In section [71 we put forward the rigorous approach relying on 
the solution of coupled integro-diffcrential equations. The differences between the two 
models are quantified in section [SI Finally, in section [31 we elaborate on the impacts 
of our work on the design and performance of frequency conversion processes for 
quantum enhanced applications. The same reasoning is then applied to the process of 
parametric down-conversion in sections[TOlonward. Section[T8lconcludes the paper and 
summarizes our findings. [Appendix A| to [Appendix F| contain additional information 
and further calculations. 

2. FC: Overview 

A general frequency conversion process is sketched in figure [U Mediated by the 
nonlinearity of the crystal and a strong pump beam two input fields a'™^ and c'™^ 
are interconverted into two output fields a^""*^ and c^°"*-'. This frequency conversion 
process is more commonly known as sum frequency generation (SFG), when the input 
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beam in combination with the pump beam generates an output field at a higher 
frequency Wout — ^in + (^p: or difference frequency generation (DFG), when a field 
with frequency ujout = '^in — is created Q 

The distinction between sum-frequency and difference frequency generation arises 
via the input wave which is fed in either the a^™' or c^'"-' port. In the scope of this 
paper a*-'"-' — >■ c^°"*^ depicts a sum-frequency generation process and c*^"'-* ^(o"*) 
labels difference frequency generation. Each crystal configuration always supports 
both processes simultaneously and we refer to the overall system as frequency 
conversion. 




out) 



Figure 1. Sketch of the frequency conversion process. Mediated by ttie strong 
pump field and tiie x'^^-nonlinearity of tiie medium parts of the fields a^'"' and 
c™ are interconverted into the two output fields a(°"*) and c(°"*). 

In this paper we go beyond the standard monochromatic single-mode description 
of frequency conversion and consider a multitude of frequencies interacting with each 
other during the frequency conversion process. This becomes especially important 
for pulsed frequency conversion experiments where wave packets, spanning several 
nanometers in width, interact with each other. In the following chapters we derive 
the properties of this transformation and compare the accuracy of different theoretical 
models. 

3. FC: Hamiltonian 

We first define the electric field operators of an optical wave inside a nonlinear medium 
as [33] 

= » ( 47rgoln3(fco) ) / ^fc ^uj{k)cxp [i{kz ~ uj{k)t)] a{k), (1) 

where A labels the transverse quantization area in the material [33] ■ We use the 
slowly varying envelope approximation, i.e. the bandwidth Alo of the considered 
electric fields is small compared to their central frequency loq (Auj <C coq) and hence 
treat the dispersion term in front of the integral n{ko) in ([1]) as a constant, using 
the value at the central wave vector kg. This approximation is justified since, in the 
remainder of this paper, we only consider electric fields not too broad in frequency, 
compared to their central frequency, and take into account the rather flat dispersion 

X In classical nonlinear optics, DFG is understood as a stimulated process. The bright pump field has 
the highest frequency and the process is seeded with a weak input field, which is enhanced through 
continuous conversion of pump photons. We, in contrast, assume a weak input field, which has the 
highest frequency and the 'seed' is the bright pump field (see | 18| for details). 
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in nonlinear crystals. Finally a{k) is the standard single-photon annihilation operator 
which destroys a single-photon in mode k and obeys 

a{k),a{k')''] ^S{k-k') [a(fc),a(fc')] = 0. (2) 

In this paper we restrict ourselves to electric fields in one dimension, this means 
we assume a fully collinear propagation of the interacting fields along one axis in a 
single-spatial mode, since a three dimensional treatment does not offer further physical 
insight into the properties of the process and clutters the calculations. 

The Hamiltonian of the frequency conversion process consists of two parts: 

(x) 

The Hamiltonian Hq ' describing the free propagation of the electromagnetic waves 
through the medium for each of the involved fields [32] 

Ht\t) = J dzEi-^{z,t)Ei+\z,t), (3) 

where the factor in front of the Hamiltonian appears due to the normalization of the 
electric fields operators in ([T]) [32]. The interaction Hamiltonian of the frequency 
conversion process is given by [35] [36] [HI [19] 

i?f ^) it) ^eo Jdz (z) (z, (z, t)Ei-Hz, t) + h.c. . (4) 

E^'^\z,t) labels the pump field driving the frequency conversion process and 
Ei'^\z,t),Ec \z,t) are the two fields which are interconverted. In the derivation 
of this Hamiltonian we used the rotating wave approximation and hence only regard 
the frequency conversion terms of the nonlinear optical process while neglecting 
the parametric down-conversion and further nonlinear interactions. The pump field 
driving the frequency conversion process is a strong optical wave. We hence treat it 
as a classical field: 

eI+'> {z,t)^ E^-^iz, t)=Ap J dka [uj{k)] exp [i{kz - Lo{k)t)] (5) 

Here Ap labels the pump amplitude and a [uj{k)\ its spectral distribution ranging from 
6{uj — ujc) for cw-laser sources up to more complicated forms in the case of pulsed laser 
systems. We further assume that the pump field is not depleted during propagation 
through the crystal since only a minor part of the strong pump beam is lost during 
the frequency conversion process. The interaction Hamiltonian ([4]) becomes: 

i?f ^) {t) ^eo [dz x(^) (z) 4+) (z, t)Ei+'> (z, t)Ei-\z, t) + h.c. . (6) 



Combining ([3]) and ([6]) the process of frequency conversion is described by the overall 
Hamiltonian: 

HMt) = Ht\t) + Hl,^\t) + Hf'^Xt) (7) 

There are a variety of different constants involved in the definition of the frequency 
conversion Hamiltonian in ([7]) (see ([I]), (O, ([5]) and ([6])). However most are dependent 
on the initial definitions of the interacting fields. In the remainder of this paper we 
merge all of them into a coupling value depicting the overall efficiency of the frequency 
conversion process rendering the presented calculations independent of individual 
notations. 
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4. FC: General properties 

The frequency conversion Hamiltonian in ([7|) generates the unitary transformation 



Ufc = Texp 



(8) 



In ^ the time-ordering operator T is crucial, because the electric field operators in 
Hpcit) a-i'C time-dependent. In turn, the Hamiltonian does not commute at different 
points in time, which renders finding a solution difficult. The structure of ^ already 
gives valuable insights into the properties of the system. The Hamiltonian in ([5]) is 
bilinear in its electric field operators — the pump is treated as a classical wave — and 
the solution hence takes the form of a linear operator transformation [37l [38l [39] 1§| In 
the monochromatic picture the solution is of the form |18j : 
a(°"*)(a;) = cos(r)a(") (w) + sin{r)c^'"\uj') 

giout) (^/) ^ cos(r)c(™) {uj') - sin(r)a(") (lu) (9) 
A single frequency a^*"' (uj) has a certain conversion amplitude sin(r) which defines 
which part of the beam will be converted to c'-°"*^(a;'), via sum-frequency or difference 
frequency generation, whereas the rest is transmitted without change to a'°"*^(cj). 
This mathematical structure is very similar to the transmissivity and reflectivity of a 
beam-splitter |40| . In the frequency conversion case, however, not two spatial modes 
but two optical modes at different frequencies are coupled with each other. 

This behavior is not immediately visible in the multimode regime, where we have 
to consider the conversion of many frequencies: 



c(°"*)(w)= / dcj'C/,(cj,cj')c(™)(w')- / dLu'V,{Lu,Lu')a^'''\Lu') (10) 



Here the functions Ua/c{'^,oj') in ([TO)) define which parts of the different frequencies 
of the input beams pass the crystal unperturbed, whereas the Va/d^T^') functions 
give the portions of the beams which are converted via sum frequency or difference 
frequency generation. 

In order to unravel the underlying structure we use the constraint that the 
frequency conversion process is a unitary transformation and hence (jlOp has to form 
a canonical transformation |38[ I39j . This imposes several conditions on the properties 
of the solution which we study in detail in [Appendix A[ Under this constraint and 
with the help of the Bloch-Messiah decomposition we rewrite (fTO]) as 

4°"*)=cos(r,)if)+sin(r,)cf) 

c'~r^=co.{r,)ct^-.Hru)Af\ (11) 
where Ak and Ck are broadband single-photon destruction operators [41] defined as: 

4°"*) = / dc.^fc(w)a(°"*)(a;) C^.""*) = / dc. ^(c.) c(°"*) (c.) 



)a(°"*Hw) 


^{out 


a(™)(w) 





4") = / dc^V^fe(c^)a(™)(w) Ci"'^^ du:M^)c'^'"H^) (12) 



§ Solving the process via the Heisenberg equation of motion (see section [7]l, yields a linear operator 
valued differential equation. Linear operator valued differential equations are solved by linear 
Bogoliubov transformation, i.e. the solution will be of the form lllOI I. Furthermore with the help 
of the Bloch-Messiah reduction, discussed in [Appendix A| we are able to rewrite the general linear 
Bogoliubov transformation | 39| into the form depicted in 111 It . 
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This means that broadband frequency conversion is equivalent to the transformation 
of an orthonormal set of optical pulses given by the broadband mode shapes in 
(|12p . Each individual pulse is, independently from the rest, partially converted and 
partially transmitted. According to ([TT|) we consequently are able to regard the overall 
frequency conversion process as a quantum pulse gate [lH [TH [20] which, depending on 
the efhciency of the process, transmits the incoming pulses unperturbed or switches 
them via frequency conversion. 

The crucial parameters of this transformation are firstly the conversion amplitudes 
rfc which give the efficiency of the process — conversion efficiency = sin^(rfc) — and 
secondly the broadband mode shapes ipk{uj), Cfc('^)) V'/c(w) and 0fc(w) that define the 
range of frequencies which are interconverted. 

5. FC: Time-ordering approximation 

In order to obtain some first insights into the frequency conversion process we use 
perturbation theory. In the case of a time-dependent Hamiltonian a Dyson series 
expansion is necessary [42] . This approach was studied in detail by Agata M. Branczyk 
in [30] and [43]. Here we shortly recapitulate the general approach, starting with the 
Dyson series expansion of ([8]): 



In our case we regard the system long after the interaction in the crystal has taken 
place and consequently use plus and minus infinity as the overall bounds for the time- 
integration. 

Unfortunately it is not trivial to perform the Dyson series expansion for FC, 
because coupled time integrals appear. To simplify this evaluation we introduce 
the time-ordering approximation. Hereby we approximate ([S]) by ignoring the time- 
ordering operator T- By dropping T in ([S]) a simple Taylor series expansion becomes 
sufficient: 



Note the differences in the integration ranges and the expansion coefficients in (|13p 
and (|14[) when terms of second order and higher are considered. 

To investigate the differences between the two approaches we explicitly calculate 
the perturbation series up to second order. For this purpose it is beneficial to work 
with the electric field operators of ([TJ in the ^-representation }33) : 




(13) 




(14) 



^(+)(z,<) = E'^-^\z,t) 
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, , du! \/ujexp[i{k{uj)z — ujt)]a{uj) (15) 
^4'!T€ocAn(uJo) J J 

We also perform our calculations in the interaction picture, this means we move into 
a new reference frame where the effects of free propagation arc not present and hence 
do not need to consider the free propagation Hamiltonians. Finally wc assume a 
crystal featuring a constant x''^^-nonlinearity extending from — ^ to ^. The frequency 
conversion Hamiltonian from ^ takes the form 

H^^''\t) = eoX^'^ [[ dzEj+Hz,t,)Ei+\z,t,)Ei~\z,t,) + h.c. 

(2)r [a [a [a I \ ■ f ^k{^^p,i^a,i^c)L 

= eox^ 'L I dup j du)a j du ca[ujp) sine I ^ 

X exp (— zAwt) a{LOa)c^ {i-^c) + h.c. 

(16) 

where we introduced the two abbreviations Afc(cL)p, Wq, u^) = kp{ujp) + ka{uJa) ~ ^c('^c) 
and Aw = ujp + cua — uJc. 

The approximations performed when ignoring time-ordering effects are located in 
the time-integrations in the term exp {—lAut). The time-integration of the first-order 
term gives identical results for both the Taylor and Dyson series expansion: 

exp (-lALot) = 27r(5(Aw) (17) 

From a physical point of view the obtained delta function 6{Auj) describes the energy 
conservation between the three interacting photons: the input, the pump, and the 
output photon. 

In the second order expansion the time integrals in two perturbation approaches 
start to deviate from each other. The time-integration of the Taylor series in ([Til) 
yields 

/oo poo 
dti / dt2 exp (— zAw^i) exp (— zAw'^2) 
oo J —oo 

= 5{Auj)S{Auj') (18) 

which is just two times the result of the first-order expansion. In the Dyson series 
however we obtain a different result [301: 



/oo nti 
dti / dt2 exp [— iAwii] exp [— iAti;'i2] 
-oo — oo 



= 2Tr^-d(Auj + Auj')(d(Auj-Auj') + -- ^——] (19) 

\ TT Aw — Auj' J 

This term shows the crucial difference between the inclusion and the exclusion of time- 
ordering effects. When we exclude time-ordering effects the time-integrations in every 
order yield delta functions (5(wp-|-cja — Wc) between three different frequencies. We can 
interpret this mathematical expression by a physical process where only three photons 
interact with each other. This means that even in the case of a multitude of photon 
conversions take place in the overall process the energy is still conserved in each photon 
triplet. Effectively this means that the spectral properties of the interacting photons 
are identical in every order of the perturbation series. 
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This is quite different from the Dyson expansion. Here delta functions of the 
form 6{Auj + Auj') appear in the second-order expansion. From a physical point of 
view this delta function describes energy conservation which does not occur between 
the three but between six photons. Either the SFG (DFG) of two photons takes place 
simultaneously or a SFG (DFG) with successive DFG (SFG) conversion of a single- 
photon takes placelil This enables an enhanced flexibility in the available energy range 
and modifies the solution, especially its spectral properties, with respect to the Taylor 
expansion. Since the energy, in these terms, only has to be conserved for the overall 
process and not the individual photon triplets the constraints imposed by energy 
conservation are relaxed. A graphical representation of this effect in the third-order 
expansion is shown in figure [2] 

(a) SFG DFG SFG 




No Time-ordering 




Figure 2. Schematic diagrams depicting tlie energy conservation during 
successive SFG and DFG processes neglecting time-ordering effects (Analytic 
model) (a) and including time-ordering (Rigorous model) (b). Neglecting time- 
ordering effects the energy is conserved between the three photons in each 
individual SFG/DFG process, whereas the rigorous model including time-ordering 
effects only requires energy conservation during the overall process [30) . 

The time-ordering approximation will consequently give accurate results as long 
as the conversion efficiency is not too high and the first-order perturbation expansion 
is the dominant term in the solution. 



6. FC: Analytic model excluding time-ordering effects 

The huge advantage of neglecting time-ordering effects in frequency conversion is the 
fact that it enables us to build an analytical model of frequency conversion, which 

II The individual orders n of the perturbation expansion do not directly correspond to the conversion 
of n photons, multiple back and forth conversions of a single-photon can appear as well. In order to 
obtain the full dynamics of the process all orders have to be considered. 
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is highly beneficial since it enables quick and straightforward access to the process 
properties: 



Ufc ^ cxp 



dtH 



(FC) 



it) 



(20) 



This formula is identical to ^ except that we dropped the time-ordering operator 
T and work in the interaction picture. It enables us to directly perform the time- 
integration in the exponential function. As in the previous section we perform our 
calculations in the cj-representation and again we assume a crystal featuring a constant 
x'-^-'-nonlinearity extending from — to ^. After a straightforward calculation we 
obtain 

. L 



dtH 



{FC) 



dteox^^^ / ^dz E^+\z,t)El+\z,t)Ei-\z,t) + h.c. 
= 5/d../d..a(..-..) 

/ Ak{uJa,UJc)L \ 

I I a[UJa)c' (LOc) + h.c. 



X smc 



(21) 



where we merged all constants into the overall factor B and Ak{LOa,uJc) ~ kp(ujc — 
Wq) + ka{i-Ua) — kc{uJc). Dctails of this calculation are given in |18| . Defining 

f[uja,i^c) ^ Ba(ujc-uja) sine I I, (22) 

the unitary transformation generated by (pO)) takes on the form: 



Ufc = exp 



dUa / dUcf{uJa,i^c)a{0Ja)c''{uic) + h.C. 



(23) 



With the help of the singular- value-decomposition theorem |44] we recast this solution 
in the broadband mode formalism presented in At first we diagonalize the 

Hamiltonian by decomposing the exponent in ((23|). via a Schmidt decomposition, as: 



^/(wa,Wc) = ^(-rfc)V'fc(a;a)(^fc(wc) 

k 



(24) 



Here both {^ki'^a)} and {(/)fe(wc)} each form a complete set of orthonormal functions 
and Tfc G T?,"*". Employing equation (j24p we rewrite the unitary frequency conversion 
transformation in (l23l) as: 



Ufc = cxp 



^(-''fc) j A^a ■>pkiuJa)a{uJa) j dUc(t>l{iOc)c^ {uJc) 



+rk / dujaik*k{^a)a^ {i^a) / du) c(l)k{^c)c{^c 



(25) 



With the help of the broadband mode operators defined in ([12)) it takes on the form 



Ufc = exp 



k 

exp 



Y^{-ru)[A,Cl-AlCu 
[{-ru) [AuCl - AlC^ 



(26) 
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The solution of ([25]) is well known and, in the Heisenberg pictures, it reads [57j : 
4°"*)=cos(r,.)if^+sin(r,)C(:") 

Cr*)=eos(r,)cf)-sin(r,)if) (27) 

This simplified analytic model features exactly the structure required by the canonical 
commutation relations discussed in sectional Only the additional fact that the input 
modes and output modes in this simplified model arc always of identical shape differs 
from the general solution pT|) . 

It is evident that this treatment ignoring time-ordering effects enables a 
straightforward analytic solution of the frequency conversion process. In contrast 
to the perturbation approach in section[5]we do not require to evaluate a perturbation 
series up to infinite order and, in some cases, even the SVD can be performed 
analytically and hence no computational effort is required at all [45]. This enables 
the efficient engineering and design of frequency conversion processes as long as the 
applied approximations hold. 



7. FC: Rigorous theory including time-ordering effects 



In order to obtain a rigorous solution of frequency conversion we have to include the 
effects of time-ordering and the corresponding relaxed energy conservation conditions 
into our calculations. For this purpose we change our analysis method and regard the 
frequency conversion process in the Heisenberg picture. This approach was already 
utilized for frequency conversion in nonlinear optical fibers in |14[ I20j and is common 
for parametric down-conversion [SlIMlllZlllHlllSlSZlllSlllSllSn]. In order to solve 
the corresponding Heisenberg equations of motion we adapt the work of Mikhail I. 
Kolobov on PDC in |32| to frequency conversion. It enables the complete time-ordered 
solution of the frequency conversion process without having to perform the Dyson 
series expansion up to infinite orders. As a first step we redefine the electric field 
operators in ([T]) according to [32] 



a{z, t) 
c{z,t) - 



1 



\/27rfco 
1 



-\/27rfco 

The Heisenberg equation of motion for a(z,t) reads 



dfc ycJ^fc) exp [i {kz — uj{k)t)] a{k) 
dfc \/ uj{k) exp [i {kz — aj(fc)t)] c(fc) 



(28) 
(29) 



dt 



a{z, t) 



HFc{t),d{z,t) 



(30) 



We first calculate the commutator of a(z, t) with the free space Hamiltonian H^°''^\t). 
The commutation relation with the free space Hamiltonian -ffp'^'' (i) directly results in 



zero whereas the commutator with Hq [t) yields [32] : 



hko / dz'a^ {z' ,t)d{z' ,t), d{z,t) 



^iko J dz' [a'^{z',t),a{z,t)]a{z',t) (31) 

The crucial part in (|3T]) is the commutator between the electric fields, at the same 
time, but at two different points in the medium: 



[a{z,t),d\z' ,t)'\ dfc ^ exp [zfc(; 



(32) 
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In order to solve the integral in p2p we Taylor expand the dispersion relation up to 
first-order and obtain: 



[a{z,t),d''{z',t)] 



27r fco 



dfc 1 



ko 



ko 



cxp [ik(z — z')] 



f dk — exp [ik{z — z')] 

2tx fco J fco 



— % 

dfc — — exp hk(z — z')] 
ko az 



27r fco 

—icd 
fco no dz 

Using ((33|) the free space commutation relation in (|3T|) evaluates to: 



-S{z-z') 



Ht\t),a{z,t) 



c d 



a{z,t) 



(33) 



(34) 



no dz 

Next we calculate the commutator of a{z, t) with the nonlinear interaction Hamiltonian 
Hj^^\t) in dH]). To simplify this problem we treat uj{k) as a constant w(fco) using the 
slowly varying amplitude approximation. The commutator consequently evaluates to 
a delta function: 

c 



[d{z,t),a\z',t)] = —6{z-z') 
no 



We obtain: 



Hf''\t), aiz,t) 



D / dz' E^-^^z' ,t)a{z' ,t)c\z' ,t) + h.c. ,a{z,t) 



D* / dz' [a^z',t),a{z,t)]El~\z',t)c{z',t) 



-D* 

h 



dz'E'f\z',t) c{z',t) (^-^^ S{z-z') 



= -l—D* E\-\z,t) c{z,t) 
nno 



(35) 



Here we merged all constants into the overall variable D. The Heisenberg equation of 
motion in pop then reads: 



^a(z,t) = - — -^a{z,t) - '-—D*E^f\z,t)c{z,t) 
at no oz nno 



(36) 



Since we are interested in the spectral distribution of the output we perform a Fourier 
transformation into the frequency domain. We regard the created electric fields long 
after the interaction and therefore integrate the time dependence from minus to plus 
infinity, as already performed in the analytic solution in section [5] 



-d{z,t) 
iijja{z, bj) 



c d 
no dz 



a(z, t) 



h no 



D*E^-'>{z,t)c{z,t) 

— —a(z,uj) — — — D* [ duj'El'^\z,Lu' — uj)c(z,u!') (37) 
no oz nno J 
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Wc further separate the interaction from the propagation effects by transforming our 
operators into the interaction picture, similar to the analytic solution in sectionlHl For 
this purpose we introduce the electric fields: 

ea{z,uj) = a(z, a;)exp [-ika{uj)z] 
ec{z, Lu) = c{z, uj) exp [—ikc{uj)z] 

e^~\z,uj) ^ E^'^\z,Lu) exp[ikp{uj)z] (38) 

With the help of ((38)) we transform the differential equation in (|37)) into the interaction 
picture: 

d if 

— €a{z,u})= - jD* duj'e^~\z,uj' - u})ec{z,uj') 

X exp {—I [kp{uj' — cj) — kc{uj') + ka{uj)] z} (39) 

In a similar manner we calculate the differential equation of the second beam and 
obtain: 

§-/c{z,Oj) = - '-D j dc.'6(+)(z,C.~^')ea(^,^') 

exp {i [kp{uj - uj') + ka{uj') - kc{uj)] z} (40) 

Finally we introduce the abbreviation Ak{u!,u!') = kp{uj' — lu) — kdui') + kaito). The 
two coupled differential equations describing the frequency conversion process read: 

-^ea{z,uj) = -^D* J duj'e^p~\z,uj' - w) exp [~iAk{iAj,uj')z] ec{.z,Lo') 

^ec{z,io) = -^D J dij'el+'>{z,uj-Lj')cxp[tAk{uj',u)z]ea{z,w') (41) 

Note the inversion of the lo and uj' factors in the Ak and functions between the two 
formulas. By defining 

/(w,w',z) = -^D*el-\z,io' -Lj)cxp[-tAk{uj,uj')z] (42) 

we may write (j4ip in a more compact notation: 

— iaiz,u})= / du}'f{u},u}',z)eciz,u}') 

^^ec{z,Lu) = - J du;'r{uj',Lu,z)ea{z,iu') (43) 



7.1. Solving the differential equations 

In order to obtain the dynamics of the frequency conversion process the differential 
equations in (|43| have to be solved. Usually operator valued differential equations 
cannot readily be evaluated and, in the case of FC, this is complicated by the fact 
that we have to solve integro-differential equations, since we consider the conversion 
of many frequencies simultaneously. Note, that (|4ip is linear in its operators and 
hence classical solution methods like the split-step Fourier inversion method have 
been applied [20l [27l [28] . In the special case of a cw-pump laser the integral in (|4T|) 
vanishes and it is even possible to find analytic solutions |32| . 

In this paper we apply a different approach — introduced by Wolfgang Mauerer 
in |46| — exploiting the fact that the structure of the solution is already known: it 
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is a linear operator transformation ^TU\\ . Using PU)) and we obtain two pairs of 
classical integro-difFcrcntial equations ^\ 

^Ua{z,L0,L0") = - j duj'f{LO,u\z)V,{z,U:',u") 

^ Vc{z,uj,u:") = f duj'r{uj',L,,z)Ua{z,uj',uj") (44) 



d 



z 



^^U,{z,u;,u;") = - / duj' f* (u;' z)Va{z,u;' ,u;") 
d 



and 



^^Va{z,uj,iu")^ j duj',f{uj,J,z)Uc{z,J,oo") (45) 

which cover the complete dynamics of the frequency conversion process. 

We solve these coupled intcgro-differential equations using an iterative approach. 
For the pair in (j44|) this means we first formally integrate both differential equations 
along z, where we assume a medium of length L, as in in the analytic solution discussed 
in section [SJ 

Uaiz.w.J') ^ 5{uj -uj') - r dz I duj'f{uj,uj\z)Vc{z,uj',uj") 



V,{z,u,uj")^ dz du'r{oj\Lo,z)Ua{z,uj',u"). (46) 



Here we also included our knowledge about the initial solution. If no interaction 
takes place the process is described by the identity operation a(°"*)(aj) = a'^™^(w) and 
g(o«t) ^ ^y^^^Yl follows: 

Ua{z,UJ,Uj") = Uc{z, UJ.Uj") = 5{uj-uj") 

Va{z,UJ,Uj")^Vc{z,UJ,J')^{) (47) 

Starting with the initial solution for Ua{z, tu, tu") we then perform the two integrations 
in (|46p and obtain a preliminary Vc(z,w,w"). This is then used to obtain a new 
Ua{z , uj , uj") . We repeat this iterative procedure till the functions converge. 

The same method is applied to the second set of differential equations defining 
Uc{z,uj,u}") and Va{z,uj,uj") which gives us the complete time-ordered solution of 
the frequency conversion process. The implementation of this algorithm is discussed 
in [Appendix D[ where we also elaborate on the numerical accuracy of the applied 
method. The program code is available on our website^ 



8. FC: Comparison betv^reen simplified analytical and rigorous approach 

In sections |6] and [7] we presented two models describing frequency conversion taking 
into account the higher-order moments: a simplified analytic model excluding time- 
ordering effects and a rigorous numerical model including time-ordering. While the 
analytic method is straightforward and only requires minimal computational effort, 
it is only an approximation to the rigorous model. In this section we compare these 
two approaches in order to investigate the validity of the simplified theory for the 



1 The program code can be downloaded from the pubhcations section on our website. The current 



url is |http : //physlkwww.unl-paderborn. de/ag/ag-sllberhorn/publlcat ions .html | 
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modelling of frequency eonversion proeesses and the effects of time-ordering present 
in the rigorous solution. 

From the discussion in section [5] we are able to directly deduce that the analytic 
solution features constant mode shapes connected to conversion amplitudes which 
oscillate with a sinusoidal pattern. In the low gain regime, when a first-order expansion 
is sufHcient, the rigorous solution of section [7] gives the identical result. However, as 
soon as the higher-order moments yield a significant contribution to the frequency 
conversion process, the rigorous solution will feature a reshaped frequency spectrum 
due to the relaxed energy conservation condition discussed in section [SJ 

Analytic model 




Figure 3. Comparison of tlie conversion amplitudes 11/0(0;, cj')] between 
the rigorous and the simphficd analytical model of a correlated frequency 
conversion process, with rising conversion efficiencies from (a) to (c). For low 
conversion efficiencies depicted in (a) both models yield identical results. When 
conversion efficiencies of unity are approached in (b) the rigorous model starts 
to deviate from the analytic approach due to the relaxed energy conservation 
condition. Significant differences occur when conversion efficiencies beyond unity 
are considered (c). 

In order to quantify the discrepancies between the two approaches we simulate two 
frequency conversion processes: a correlated frequency conversion where many optical 
modes are converted and an almost uncorrelated frequency conversion process where 
only a few optical modes rk are present |18| . The frequency conversion process is 
defined by the incoming pump laser and the applied nonlinear crystal, i.e. the length 
of the material and its dispersion properties. These parameters can be represented 
by just six values: The length of the crystal, the group velocities of the three beams 
and the width and the amplitude of the Gaussian pump beam of our ultrafast laser 
system driving the frequency conversion process. The details of this procedure and 
the applied process properties are given in [Appendix C[ 

The crucial parameter in the solution are the '^O | amplitudes in the 
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Analytic model 




Figure 4. Comparison of the conversion amplitudes | V'a(i^, i^') | between 
the rigorous and the simplified analytical model of an uncorrelated frequency 
conversion process, with rising conversion efficiencies from (a) to (c). For low 
conversion efficiencies in depicted in (a) both models yield identical results. When 
conversion efficiencies of unity are approached in (b) the rigorous model starts 
to deviate from the analytic approach due to the relaxed energy conservation 
condition. Significant differences occur when conversion efficiencies beyond unity 
are considered (c). 

general solution in (1101) which describe the conversion amplitude — conversion 
efficiency = jconvcrsion amplitude^ — at individual frequencies. The 
amplitudes for the correlated frequency conversion case are presented in figure [3l 
Figure [Hla) gives |ya(a;, a;')| at a low conversion rate. In this regime both models 
yield identical results. At efficiencies approaching unity, as depicted in figure [3] (b), 
minor deviations between the two models appear. Figure [3] (c) shows |ya(w,a;')| for 
a conversion efficiency beyond unity. In this regime the time-ordering approximation 
breaks down and the two theories predict significantly different results. While the 
analytic model predicts a dip in the conversion efficiency about the central frequencies, 
the rigorous model still shows a Gaussian shape. 

For the case of uncorrelated frequency conversion the |VQ(a;,a;')| amplitudes are 
given in figure U) As in the correlated case|4](a) presents the process at low conversion 
efficiencies where both models are identical, (b) depicts near unit conversion efficiency 
with minor deviations between the two approaches and (c) visualizes the discrepancies 
between the two models arising at conversion rates beyond unity. In this uncorrelated 
frequency conversion case the discrepancies between the different models are much 
more prominent than in the correlated frequency conversion case. In the analytic 
model the side peaks are much more prominent as in the rigorous model. Furthermore 
the rigorous model shows a broadening and along the +45° axis, which is a direct result 
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of the relaxed energy conservation conditions in this model. 

The corresponding frequency conversion amplitudes sin(rfc) and mode shapes 
are displayed in figure [S] for correlated and in figure IHl for uncorrelated frequency 
conversion. The figures in the column on the left show the conversion amplitudes 
sin(rfc) — conversion efficiency — sin^(rfc) — whereas the two columns on the right 
present the corresponding mode functions ipi{iy),ipi{i>), and ^i(i^) for the first 

optical mode featuring the highest conversion efficiency, where and tpi{^) 

well as '^i(i') and ^i(i') arc of identical shape. 
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Figure 5. Comparison between tlie rigorous and the analytical approach in 
correlated many mode ultrafast frequency conversion. In the low conversion 
regime, presented in (a) (15%/15% conversion efficiency in the first mode in 
the analytic / rigorous model), both approaches evaluate to identical results. 
Approaching unit efficiency in (b) (86%/86% conversion efficiency in the first 
mode) the two approaches start to show differences, which become significant 
when optical gains beyond unity arc considered in (c). (100%/99% conversion 
efficiency in the first mode) 



In the conversion regime up to unity (a) both approaches yield almost identical 
conversion amplitudes sin(rfe) and corresponding mode functions ipi{i'),ipi{i'), (pi^v) 
and ^lif)- At conversion rates about unity minor discrepancies in the mode shapes and 
sin(rfc) amplitudes start to appear. Significant differences between the two theories, 
however, are only present when we consider rates beyond unit conversion efficiency as 
depicted in (c). 

Interestingly the time-ordering effects affect the correlated and uncorrelated 
frequency conversion in a different manner: in the correlated case the rigorous 
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Figure 6. Comparison between the rigorous and the analytical approach in 
uncorrelated few mode ultrafast frequency conversion. In the low conversion 
regime, presented in (a), both approaches evaluate to identical results (6.4%/6.3% 
conversion efficiency in the first mode in the analytic / rigorous model). 
Approaching unit efficiency in (b) the two approaches start to show differences 
(100%/89% conversion efficiency in the first mode), which become significant when 
optical gains beyond unity are considered in (c) (30%/99% conversion efficiency 
in the first mode). 



model predicts a considerably faster drop of the conversion amplitudes sin(rfc) than 
the analytical approach. Diamentral to this behaviour the sin{rk) amplitudes of 
decorrclated case in the rigorous model, presented in[Bl^c), remain at unit conversion 
efficiency once they reach this value. In general it is not possible to associate the time- 
ordering effects with a specific behaviour of the broadband conversion amplitudes 
sin(rfc). Consequently each frequency conversion process has to be investigated 
individually. 

Both correlated and uncorrelated down-conversion, however, share the property 
that the rigorous model predicts a significant broadening in the corresponding mode 
shapes in the high gain regime. Please note, that in figure [5] (c) the first mode of the 
analytic model is not the Gaussian from [5] (a) and (b) any more since we always plot 
the modes with the highest conversion efficiency and, in this power regime, this is not 
a Gaussian any more. 

In summary the analytical model accurately describes frequency conversion in 
the low gain regime up to unit conversion efficiency, where minor deviations start 
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to appear. Only in the high gain regime, when conversion efBcieneies beyond unity 
are reached, comphcated non-trivial deviations from the analytical model have to 
be taken into account and a rigorous treatment of frequency conversion is required. 
For most experimental setups and applications it is hence perfectly justified to apply 
the simplified analytic model to minimize the computational effort, as long as its 
limitations are kept in mind. 

9. FC: Outlook 

These time-ordering effects have practical implications for the deployment of frequency 
conversion processes. Especially affected are frequency conversion processes which 
serve as quantum pulse gates [THl [HI HO] • In theory a perfect quantum pulse gate 
converts a single optical mode with unit efficiency. Actually building a quantum pulse 
gate, however, is complicated by time-ordering effects. As is evident from figure [H] the 
relaxed energy conditions move the frequency process from the single-mode regime 
towards a more multimode behaviour. This effect fundamentally limits the gate 
performance at high conversion efficiencies. One could, in principle, use advanced 
engineering techniques, such as hypergrating structures to reduce the amount of 
multimodeness in the state [51j . yet still the time-ordering effects seem to remain a 
fundamental limitation. Whether or not it is actually possible to engineer single-mode 
pulse gates including the effects of time-ordering remains an open research question. 

10. PDC: Overview 

Having discussed frequency conversion we now extend our investigations to the second 
important nonlinear optical process: parametric down-conversion. The principle of 
PDC is sketched in figure [T] A strong pump beam decays inside a nonlinear optical 
material into two beams commonly labelled signal and idler. The output state formed 
by the two beams is known as finitely squeezed EPR-state or twin-beam stat43- As 




Figure 7. Sketch of the parametric down-conversion process. A strong pump 
beam decays inside the nonhnear optical material into two beams usually labelled 
signal a and idler b, forming a finitely squeezed EPR-state. 

in the frequency conversion case, we do not restrict ourselves to a discussion of the 
parametric down-conversion process in the monochromatic picture, but extend the 
theories to the multimode picture including the interaction of many frequencies. This 
is especially important when the PDC process is pumped by pulsed laser systems. 

+ PDC is also able to create squeezed states. This process is discussed in |27ll28j . 
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The interaction Hamiltonian of the PDC process — using electric field operators as 
defined in ^ and the rotating wave approximation — takes on the form: 

7?f^^)(t)=eo j dzx'^^\z)&+\z,t)Ei-\z,t)E\r\z,t) + h.c.. (48) 

As in the frequency conversion process we assume a strong pump field exceeding the 
amplitudes of the signal and idler fields by several orders of magnitude and hence 
treat it as a classical field propagating undeplcted through the medium (see (O). This 
transforms the interaction Hamiltonian to: 

7?f^^)(0=eo j Azx^^\z)E\,+\z,t)El\z,t)E[-\z,t) + h.c.. (49) 

Using the free space propagation Hamiltonian from ([3]) and the interaction 
Hamiltonian from (|49p the process of parametric downconversion is given by the overall 
Hamiltonian: 

HpDc{t) = H^^'^Ht) + H^^\t) + ^r^^^(t) (50) 

While, at first glance, the process of parametric down-conversion seems very different 
from the process of frequency conversion, comparing the interaction Hamiltonian of 
PDC in and FC in ^ reveals that they both feature bilinear Hamitonians — the 
pump is treated as a classical field — with an almost identical structure and hence 
share many mathematical properties. 

It is therefore straightforward to extend our presented calculations corresponding 
frequency conversion to parametric down-conversion. In order to avoid repetition we 
are going to only state the results and elaborate on the differences and similarities 
to the process of frequency conversion. The detailed calculations will be presented in 



12. PDC: General properties 

The general solution of (jSOp takes on the form of a linear operator transformation 



6(°"*)(w)= / dw'L/6(w,w')^*"'(^')+ / dLj'Vb(Lj,Lj')a(")t(^') (51) 



This solution is constrained by the fact that it has to form a canonical transformation 
[3S1 [Sni HZ] • Under this restriction we are able to rewrite it as: 

4°"*' = cosh(rfc)4™) -Hsinh(rfc)Bi™)^ 

= cosh(rfe)^[:") + sinh(rfe)4'"^^ (52) 

where Aj. and Bf^ are defined as broadband single-photon destruction operators [41] : 

= Jdcuipkico) a(°"*) (lo) ^Jdu a (to) (w) 

4™^ = / dL.Vfc(^)a<"'(^) = I dc^0fc(c^)6(™)(c^) (53) 
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The details of this procedure are given in [Appendix B[ 

According to (j52[) the parametric downconversion process is a broadband twin- 
beam squeezer transformation. However, not a single finitely squeezed EPR state, 
but a multitude of twin beam states are emitted into broadband ultrafast pulse 
modes Ak and Bk- The crucial parameters of this transformation are firstly the 
amplitudes rk which give both the amount of generated EPR entanglement — 
squeezing[dB] = — lOlogj^Q (e^'^'"'') — and the number of emitted EPR states, and 
secondly the mode shapes (/3fc(a;), ^fc(i^); ^/'fe(w), and (pkX^) which define the form of 
the optical modes in which the states are created. 



13. PDC: Time-ordering approximation 

To obtain some first insights into the parametric down-conversion process given by 
(|50p we applied perturbation theory, as performed in the case of frequency conversion, 
discussed in section [5] The mathematics are identical to the ones presented in section 
[5] Including time-ordering and performing a Dyson series expansion is involved yet 
gives an enhanced flexibility in the energy conservation condition, whereas neglecting 
time-ordering enables a straightforward Taylor expansion. An extensive discussion of 
Dyson series vs. Taylor series for PDC was performed by Agata M. Branczyk in [30] 
and gg. 



14. PDC: Analytic model excluding time-ordering eflfects 



As in the frequency conversion case, presented in section [BJ we first solve the process 
excluding time-ordering effects. Again we use the electric fields in the frequency 
domain (|15|) and move into the interaction picture: 



-- / dtH 



Up DC = exp 



Retracting the steps from section |6] we obtain 



(PDC) 



it) 



Up DC = cxp 



duJn 



where f{u!a^LOb) is defined as 

f{uja,^b) = B a{u!a + ^b) siuc 



and Ak{uJa,ijJb) = kp{uja + UJb) - ka{uja) - kb[uJb)- 

In the broadband mode formalism this process evaluates to 



Up DC = exp rk (ilSfc - ifcSj 

k 

= (g)cxp [rk [AlBl-AkBk 

k " k 

which, in the Heisenberg formalism ([57]) . takes on the form: 



^Sfhrk) 



(out) 



cosh{rk)Al^"^ + smh{rk)B 



t(m) 



(out) 1 / \ -ft(in) , • 1 / \ (in) 

l ' ^ cosh{rk)B}. ' + smh{rk)Ai: ' 



(54) 



(55) 



(56) 



(57) 



(58) 
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The details of this calculation are given in [T^] . 

This result exhibits exactly the structure imposed by the canonical commutation 
relation in ([5^ . except for the fact that, as in the frequency conversion case, the input 
and output modes are of identical shape. 

In conclusion the analytic model ignoring time-ordering effects enables a 
straightforward solution of the parametric down-conversion process, which enables 
the efficient engineering and design of frequency conversion processes as long as the 
applied approximations hold. 



15. PDC: Rigorous theory including time-ordering eflfects 

Having elaborated on solving parametric down-conversion neglecting time-ordering 
effects we now explicitly include these effects into our analysis. For this purpose we 
adapt the approach presented in section [71 We first define new electric field operators 
according to [52] 

1 



a{z,t) = —j= / dfc ■\/a;(fc) cxp [i (kz — uj(k)t)] a{k) 
V27rfco 

1 



b{z, t) = / dfc Vw(fc)exp [i {kz - uj{k)t)] b{k). (59) 

V27rfco J 

In terms of ([S^]) the Hcisenbcrg equation of motion for a{z,t) takes on the form: 

j^a{z, t) = j^ [HpDc{t), aiz, t)] (60) 

Repeating exactly the same steps as in section[7|we obtain two operator valued integro- 
differential equations describing the downconversion process 

d if 

—ea{z,Uj) = --D / du' cxp[tAk{uj,uj')z]E^'^\z,u + oj')el{z,u') 

—h{z,uj)^-^D dLu'ciip[iAk{u',uj)z]El+Hz,Lu + Lj')el{z,uj') (61) 

where we introduced the shorthand Afc(a;,w') = kp{u! + uj') — fca(w) — fcb(cj'). The 
structure of this result is very similar to the equations derived by [IHl HSl [IZ] which 
serves as a nice cross check of our calculations. Also please take note of the switch of 
oj and uj' in the second equation. 
By defining 

/(cj,cj',z) = -^DEI+\z,uj + tj')eyii>[tAk{uj,uj')z] (62) 

we rewrite the differential equations in a compact form: 
d f 

— iaiz^uj) = / duj'f{uj,uj',z)el{z,uj') 
d f 

— ib{z,uj)= / duj'f{uj',uj,z)e'l{z,uj') (63) 



15.1. Solving the differential equations 

Since the structure of the two differential equations describing the PDC process in 
(|63p is identical to the ones describing the frequency conversion process (|^T|) and 
we apply the same solution method as presented in section [73] 
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We obtain four classical integro-differential equations. Two for Ua{z,uj,uj') and 
Vb{z,uj,uj'): 

^Ua{z,L0,L0") - j Ajf{uj,Uj',z)V:{z,u',u") 

^ Vb{z,uj,u:") = I duj'f{uj\uj,z)U:{z,uj',u") (64) 



dz 

And two for Ub{z,uj,uj') and Va{z,LO,Lo' 
d 



g^Ub{z,uj,uj") = J duj'f{uj',u,z)V:{z,uj',uj") 

^Vaiz,io,io")^ j du'f{u:,u:',z)U;{z,oj',u;") (65) 
The initial conditions are: 

Ua{z, LU, Uj") ~ Ub{LL>, Ll>" , z) ^ S{u! — Lo") 

Va{z,UJ,Uj")^Vb{uj,Uj'\z)^{) (66) 

As in the frequency conversion case they can be solved via an iterative approach. 
Details of this calculation and the numerical errors in the solution method arc give in 
[Appendix F[ The program code is available, together with the frequency conversion 
code, on our website. 



16. PDC: Comparison between simplified analytical and rigorous 
approach 

In sections[T4]and[l5]we presented two models describing parametric down-conversion: 
a simple analytic model excluding time-ordering effects and a complex numerical 
model including time-ordering. As in the frequency conversion case the analytic 
method is straightforward and only requires minimal computational effort while the 
rigorous approach demands the solution of coupled integro-differential equations. In 
this section we compare these two approaches in order to test the validity of our 
simplified theory for the modelling of parametric down-conversion and the effects of 
time-ordering present in the rigorous approach. 

As in the frequency conversion case, presented in section |51 we consider a 
correlated and an almost uncorrelated process pumped by ultrafast pump lasers (see 
[Appendix E| for the process properties and [Appendix F for details on the numerical 
implementation) . 

The crucial parameters in the solution are the |Va/b(w,cj')| amplitudes in the 
general solution in (jSip which describe the amount of EPR-squeezing that is generated 
during the nonlinear interaction. 

The resulting |Va(aj,a;')| for a correlated parametric down-conversion process is 
depicted in figure HI The results for uncorrelated parametric down-conversion are 
visualized in figure |9l Similar to the frequency conversion case we choose to depict 
the low down-conversion regime in (a), the regime where first discrepancies occur in 
(b) and the efi^ects in the high down-conversion regime in (c) . 

The corresponding EPR squeezing values and mode shapes are depicted in figures 
1101 and [11] for rising down-conversion rates from (a) to (c). The figures in the 
column on the left shows the squeezing squeezing values obtained via the relation: 
squeezing[dB] = — lOlog^Q (e^^'''=) and the corresponding mean-photon number of the 
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Analytic model 





Figure 8. Comparison of the process amplitudes IV'aC'i', t^')! between the rigorous 
and the simplified analytical model for a correlated parametric down- conversion 
process. In the low squeezing regime in (a) both models yield identical results 
((n) = 0.59/0.60 in the analytic / rigours model.). Only at EPR squeezing 
values of 12 dB, presented in (b), minor discrepancies start to appear ((n) = 
44.14/53.18). Significant differences are only present for extremely high down- 
conversion rates visualized in (c) ((n) = 1708.16/3328.27). 



complete state, whereas the two columns on the right present the corresponding mode 
functions ipi{u),'ijji{i'), and ^i(i^) for the first optical mode, where </'i(^) ^-^d 
as well as (t)i{v) and S,i{v) arc of identical shape. "Analytic model" labels the 
solution excluding time-ordering effects, as presented in section [Ml and the "Rigorous 
model" label marks the rigorous solution from section [TSl 

Up to squeezing values of 12dB presented in (b) the two approaches give identical 
results, only when squeezing values beyond this bound are considered significant 
differences between the two models start to appear. The rigorous model predicts 
more squeezing than the analytic model and the relaxed energy conservation condition 
induces a broadening of the spectral distributions |T4/f,(w, w')| which leads to a 
broadening of the mode shapes in the high gain regime. In general the impacts of 
time-ordering on correlated PDC are much weaker than the changes present in the 
uncorrelated PDC case. 

In summary the analytical model accurately describes parametric down- 
conversion in the low gain regime up to two- mode squeezing values of 12dB, where 
minor deviations start to appear. Only for extremely high squeezing values, in the 
range of 20dB and higher, complicated non-trivial deviations from the analytical 
appear and give significant contribution which require a rigorous treatment of 
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Figure 9. Comparison of the process amplitudes | Vti(ct^, u)')\ between the rigorous 
and the simplified analytical model for an almost uncorrelated parametric down- 
conversion process. In the low squeezing regime in (a) both models yield identical 
results ((n) = 0.07/0.07 in the analytic / rigours model.). Only at EPR 
squeezing values of 12 dB, presented in (b), minor discrepancies start to appear 
((n) = 2.80/4.08). Significant differences are only present for extremely high 
down- conversion rates visualized in (c) ((n) = 39.39/279.87). 



parametric down-conversion. For most experimental setups and applications it is 
hence perfectly justified to apply the simplified analytic model to minimize the 
computational effort, as long as its limitations are kept in mind. 

17. PDC: Outlook 

The time-ordering effects on PDC have practical implementations for the deployment 
of PDC sources in quantum optical experiments. In contrast to FC, however, the 
time-ordering effects are beneficial to the performance of the sources. They lead to 
much higher EPR squeezing values as predicted by the simplified analytic model. 

18. Conclusion 

In conclusion we developed two models for the nonlinear optical processes of frequency 
conversion and parametric down-conversion taking into account higher-order photon 
number effects. The presented rigorous numerical model relies on the solution of 
coupled differential equations, whereas ignoring time-ordering effects enabled us to 
construct an analytical solution. 

Our analysis revealed that the analytic model gives accurate results for frequency 
conversion processes up to unit conversion efficiency and parametric down-conversion 
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V (a.ii..) 



Figure 10. Comparison between tlie rigorous and the analytical approach in 
correlated many mode parametric down-conversion. For low down-conversion 
rates, presented in (a), both approaches evaluate to identical results. Only in the 
case of rising squeezing values in (b), with squeezing values about 12dB, the two 
approaches start to show minor differences, which become more prominent when 
even higher squeezing values are considered (c). 



process up to 12dB of squeezing. However as soon as we go beyond these bounds, the 
relaxed energy conservation conditions, induced by the time-ordering effects, require 
a rigorous treatment of these processes. 

The developed fundamental understanding of frequency conversion and 
parametric down-conversion gives valuable insights into the quantum properties of 
the nonlinear optical processes, extends the current theoretical descriptions of FC 
and PDC, and simplifies considerably the engineering process for future quantum 
information applications. 
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Figure 11. Comparison between the rigorous and the analytical approach in 
(almost) uncorrelated few mode parametric down-conversion. For low down- 
conversion rates, presented in (a), both approaches evaluate to identical results. 
Only in the case of rising squeezing values in (b), with squeezing values about 
12dB, the two approaches start to show minor differences, which become more 
prominent when even higher squeezing values are considered (c). 



Appendix A. FC: Canonical transformation conditions 



The frequency conversion process in (jlOp is a unitary transformation, hence the 
generated output modes must preserve the canonical commutation relations. This 
imposes several restrictions on the structure of the solution. We evaluate these by 
extending the calculations from [351 131] to frequency conversion. At first we rewrite 
(jlOP in the more compact notation 

V (out) (a) - (m) (a) -.(m) 



a} 



^{out) 



(a) ^ (m) 
(c) -(in) 



(a).(i' 
(c) ^ (in) 



where i and j label the individual frequencies of the electric fields and summation over 
repeated indices is understood. These two input-output relations must preserve: 

(A.2) 









Oi, a] 




Ci,c] 







Using (jA.l[) and (jA.2[) we obtain three conditions for frequency conversion: 

UaUl + VaV} = UJJI + ^ t 



(A.3) 



Theory of quantum FC and PDC in the high gain regime 27 
UaV^ ~ VaUl = (A.4) 

Furthermore the commutation relations have to be preserved for the inverse 
transformation as weh: 

cf = ufh^;^'^ + vf^a^;^'^ . (A.5) 

With (IA.5I) the canonical commutation conditions in (jA.2[) yield the restrictions: 

UlUa + = UlUc + V^Va = 1 (A.6) 

UlVa - K^C/c - (A.7) 

The equations (|A.3p , (|A.4p , (jA.6|) and (|A.7[) impose several constraints on the solution. 
However they are rather unintuitive representations of the symmetries governing the 
frequency conversion process, yet with the help of the Bloch-Messiah reduction [31] 
it is possible to unravel their underlying structure: As a first step we decompose the 
four matrices Ua,Va,Uc,Vc as 

= A'^D^Bf Vc = (A.8) 

where A and B are unitary matrices and D is a diagonal matrix with real entries. 
This definition is equivalent to a singular- value decomposition except for the fact that 
we allow the individual elements in D to exhibit negative values^ 

The matrices UaU^ and VaV^ are Hermitian and (|A.3p implies that they commute 
hence both are diagonalised by the same unitary matrix P: 

P UaUlP^ = D P VaV^P^ = D' (A.9) 

With the help of the decomposition in ()A.8P they take on the form: 

PAlDfAl'p^ = D PA^DfAl'p^ = D' (A.IO) 

And we obtain A^ = A^. From (jA.3[) using UcU^ and VcV^ we infer in a similar 
manner A" = A". Evaluating the conditions for the inverse transformation using 
(|A.6p yields = _B" and ~ B"^ . Consequently the decomposition in (|A.8p can be 
written as: 

Ua = AaDlBl = A^D^.Bl 

Va = AaDlBl V, = A^DlBl (A.ll) 

Using the matrices in (jA.llj) in conjunction with the conditions in ()A.3P we further 
obtain: 

Dl-+Df=l Df+Uf^l. (A.12) 

This implies that the individual elements of the D matrices have to obey cos(rfc) and 
sin(ri.) behaviour. Applying the conditions in (jA.6p to the transformation matrices in 
(|A.11|) resuhs in: 

Df + Df =1 Df+ Df ^ 1 (A.13) 

2 2 2 2 2 2 

From which we conclude D'^ = = £>" and = Z?^' = D^' . Taking everything 
into account the final decomposed frequency conversion matrices read: 

Ua = AaD^Bi = A^D^Bl 

Va = AaD-Bl Vc - A^D^'Bl 

L>"' (A.14) 

* The reason for this extension of the singular- value decomposition becomes clear in I IA.I6I 1 . 
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In the original representation they take on the form: 

Ua{uj,uj') = (fil juj) COs(rfc)-0fc {uj') 
k 

Va{uj,uj') = ^Iji^) sin(rfc)(/)fc(iLj^) 

k 

Uc{uj,uj') = ^^^(a;)cos(rfe)(/)fe(a;') 

k 

V,ioj,Lj') = J2Ck{io)sm{rk)M'^') (A.15) 

k 

From these symmetry relations the frequency conversion process in pop must take on 
the form of multiple beam-splitter relations in orthogonal modes 

4°"*)=cos(r..)i^)+sin(r.)(7f^' 

C^""*) = cos(r,)cf' - sin(r,)if ) (A.16) 
where we defined: 

4"^ = J dc^V^fe(c^)a(*")(w) C^") = J dw</)fc(L.)c(™)(w) (A.17) 

Note however that the canonical commutation relations do not demand that the input 
and output modes are of identical shape. In principle the input modes A^"^^ and 
output modes C*'*"' could feature completely different spectral mode functions <^k{oj) 
and ipk{'^) but still form a canonical and hence unitary solution. 

Appendix B. PDC: Canonical transformation conditions 

As in the case of frequency conversion the parametric down-conversion process is a 
unitary transformation, hence it must preserve the canonical commutation relations. 
Retracting the calculation in ( [Appendix A[ ) for parametric down-conversion they read: 



UaUl - VaVl = u^ul - v^yl = / (B.l) 

VaV^ - VaUl = (B.2) 

For the inverse transformation they evaluate to: 

UlUa - {Vlv^f = ulUt - {VlVaf = I (B.3) 

UlVa ~ {Ulv.f - (B.4) 



With the help of the singular- value-decomposition theorem and (jB.ll) , (|B.2[) and (jB.4p 
the four matrices of a general parametric down-conversion process in (jSip take on the 
form: 

Ua{u:,u)') = ^(^^(a;)cosh(rfc)'0fc(a;') 

k 

Va{^,^') = ^(^fe(t^)sinh(rfc)0fc(a;') 

k 

Ub{uj,uj') = ^^fc(w)cosh(rfc)0fe(a;') 

k 

Vi,iu;,u;') = Y.ek{^)smHr,)rkito') (B.5) 

k 
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From these symmetry relations the parametric down-process (|5ip consists of multiple 
twin-beam squeezers in orthogonal optical modesQ 

= cosh(rfc)if ^ +sinh(rfc)Bf 

= cosh(r,)sf + sinh(r,)if (B.6) 

where we defined: 

4°"*^ = J dc.v5fe(c.)a(°"*)(w) = J dc.a-(w)S(°"*)(w) 

Note however that, as in the frequency conversion case, the canonical commutation 
relations do not demand that the input and output modes are of identical shape. 



Appendix C. FC: Simulated frequency conversion processes 

In our simulation of frequency conversion we did not restrict ourselves to a specific 
crystal material or wavelength range, but created a generic model of the process. For 
this purpose we first moved from the (w, aj')-system to the parameter range 
relative to the central frequencies of the frequency conversion process {ujq,lo'q). In the 
simulation we work with a Gaussian pump distribution, as created by pulsed laser 
systems. The pump distribution in ([2^ and (PT|) takes on the form 

a{h' — v') = Ep exp 

where Ep labels the pump amplitude and a the pump width. The second function we 
have to adapt is the phasematching function Afc(w, to') = kp{uj' — w) — fcc(w') -I- ka{ui). 
As a first step we perform a Taylor expansion of the individual fc(w) terms up to first 
order about their central frequency ujt^: 

k{uj) ^ k{ujn) + -^k{ujQ) {u} - ujo) (C.2) 

This is justified since we restrict ourselves to nonlinear processes not to broad in 
frequency (slowly varying envelope approximation Aoj <C i^o) fa-r off any singularities in 
the dispersion relation. At the central frequencies the process, per definition, displays 
perfect phasematching kp{uj'Q ~ loq) — fcc('^o) + fca('^o) = and the phasematching 
function simplifies to: 

Ak{v, v') = -^kpi^'o - ^o)i'^' -v)- ^fcc(wo) v' + ^^^(wo) v (C.3) 

The three remaining parameters ^A:p(a;g — Wo)i '^^ci^o) and -£jka{uJo) — the inverse 
group velocities of the three interacting beams — define the material properties of the 
system and can be adjusted accordingly. 

This compact notation enables us to simulate any frequency conversion process 
with the help of just 6 parameters: The width and amplitude of the pump beam, the 
group velocities of the three interacting waves and the length of the nonlinear medium. 

tl In principle the twin-beam squeezer has a phase degree of freedom which we absorb in the 
definition of the electric field operators Af^ and Bf^. 



2a2 



(C.l) 
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In order to evaluate a correlated FC with many modes r^, as presented in figure 
[3] and [51 we use the parameters a = 4.0, ^fcp('^o ^ "^o) = 0.5, -£jkc{uj'Q) = 3.0, 
^/Ca(cjo) =4.5 and a crystal of length L ^ 2. The pump amplitude Ep is adjusted to 
give the targeted down-conversion rates. 

For the (almost) uncorrelated case with few optical modes , as depicted in figure 
SandEl we applied a = 0.98190, ^fc^K-wo) = 3.0, ^k^iuj'o) = 1.5, £kaiujQ) = 4.5 
and a crystal of length L — 2. Again the pump amplitude Ep is adjusted to give the 
desired down-conversion rates. 



Appendix D. FC: Numerical implementation 

In order to obtain the time-ordered solutions we solved the classical differential 
equations in (j^J) and ps)) which give the functions Ua{z , oj , uj") , Udz , uj , uj") , 
Va{z,Lu,uj"), and Vc{z,uj,uj") describing the frequency conversion process. 

In the numerical implementation of frequency conversion we used a sampling 
of 500 points for each frequency degree of freedom and 500 points in z-direction 
to discretize the functions [/^(z, w, w"), Uc{z,ui,uj"), 14(2;, w, o;"), Vc{z,uj,uj") and 
/(oj, oj', z). We evaluated the successive integrations in pS)) via the trapezoid rule until 
the solutions converged. The actual solution defining the overall process properties is 
given by the matrices at the end of the crystal Ua{z = ■|,a;,a;"), Uc{z — ^,uj,uj"), 
K(^= and y,(z= 

We checked the accuracy of the result in a variety of ways: At first we evaluated 
the canonical transformation conditions in (|A.3[) . (IA.4I) . (jA.6[) . and (|A.7|) . For example 
in the case of (jA.4p we calculated: 

dw'K(z = ^,u,J)U,iz = ^,Lo",u'r 

= (D.I) 

and determined the distance of D^'^^^^^{z = ^,uj,uj') from the expected zero matrix 
and consequently the error in the solution via: 

fdcj rdcj'i?(*//)(z = ^,w,w') , , 

error = - - - — ■ - fD 2) 

0.5[JdujJduj'\Vaiz = ^,u;,u;')\+Jdu;Jdu;'\Uciz = ^,uj,uj')\] 

In all presented cases the obtained error was below 0.0085. 

We also checked the numerical Schmidt decompositions of Ua{z = ■j,a;,a;"), 
Uc{z = ^,uj,uj"),Va{z = ^,uj,uj"), and Vc{z = ^,u}, uj") to verify the mode properties 
derived in [Appendix A[ During this process the decompositions of Ua{z = ^^uj^uj"), 
Uc{z = ^,ijj,Lij") showed numerical issues, these however could be resolved by 
decomposing Ua{z = ^,uj,uj")W^{z = ^,uj,uj"), WJz = i^,uj,uj")Ua{z = f ,w,w"), 
Ua{z = ^,u:,u:")Ul{z = f,w,w"), C/t(z = f,cj,w")?/c(z = f,^,'^") instead. 
The obtained modes from these four matrices provided a much improved stability 
especially in the high gain regime. Using these Schmidt modes we verified that the 
obtained Schmidt values of the U and V matrices behaved like cos(rfc)^ -l-sin(rfe)^ = 1 
with errors below 0.0063. We also asserted that the decompositions yielded the 
functions (pk,'4'k,4>k-:£,k with symmetries as detailed in (|A.15[) . which were fulfilled 
within numerical accuracy. 



Theory of quantum FC and PDC in the high gain regime 



31 



The program code, published on our website, is able to direetly create the 
investigated frequency conversion processes and also performs all mentioned tests. 



Appendix E. PDC: Simulated down-conversion process 

As in the simulation of frequency conversion processes in |Appendix C| wc didn't restrict 
ourselves to a specific crystal material and wavelength range but created a generic 
model of the process. Again we first move from the (w, a;')-system to the parameter 
range {v,v') relative to the central frequencies of the parametric down-conversion 
process (ujq,uj'i^). As in the frequency conversion case we used a Gaussian pump 
distribution for the simulation, which in (1561) and (1611) take on the form 



a{v + I/') = Ep exp 



(E.l) 



2a2 

where Ep labels the pump amplitude and a the pump width. The second function we 
have to adapt is the phase-matching function Afc(cj, uj') = kp{io + io') — ka{i-u) — ki,{uj'). 
As a first step we perform a Taylor expansion of the individual k{uj) terms up to first 
order about their central frequency loq: 

This is justified since we restrict ourselves to nonlinear processes not to broad in 
frequency (slowly varying envelope approximation Aw ^ wq) far off any singularities 
in the dispersion. At the central frequencies the process, per definition, displays perfect 
phasematching kpiio^ + wq) — fc(j(wo) — fct,(wQ) = and the phasematching function 
simplifies to: 

Afc(z^, v) = ^fcp(^o + ^o)(i^' -^v)- ^fca('^o) V - -^^bii^'o) ^' (E-3) 

The three remaining parameters ^A:p(a;o -I-wq), ^fca(a;o) and -£^kh{ujQ) — the inverse 
group velocities of the three interacting beams — define the material properties of the 
system and can be adjusted accordingly. 

This compact notation enables us to simulate any parametric down-conversion 
conversion process with the help of just 6 parameters: The width and amplitude of 
the pump beam, the group velocities of the three interacting waves and the length of 
the nonlinear medium. 

In order to evaluate a correlated PDC with many modes r/j, as presented in 
figure [8l and [TOl we use the parameters a — 4.0, -^kp{ijjQ +a;o) ~ 0.5, ^fca(a;o) = 3.0, 
^A:b(wg) = 4.5 and a crystal of length L = 2. The pump amplitude Ep is adjusted to 
give the targeted conversion efficiencies. 

For the (almost) uncorrelated case with few optical modes r^, as depicted in 
figures m and HH we applied a = 0.96231155, ^kp{LJo + uJq) ^ 3.0, ^^^(wo) = 4.5, 
-^ki,{uj'Q) = 1.5 and a crystal of length L = 2. Again the pump amplitude Ep is 
adjusted to give the desired conversion efficiencies. 



Appendix F. PDC: Numerical implementation 

In order to obtain the time-ordered solutions we solved the classical differential 
equations in ([M]) and (p5|) which give the functions Ua{z,uj,uj"), [/^(z, w, w"), 
Va{z,ui,ui")^ and Vb{z,uj,uj") describing the parametric down-conversion process. 
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In the numerical implementation of parametric down-conversion we used a 
sampling of 500 points for each frequency degree of freedom and 500 points in z- 
direction to discretize the functions C/a(2;, w, w"), J7b(z, w, w"), Va{z,uj,uj"), Vf,(z, w, w") 
and f [bj , Lo' , z) . As in the frequency conversion case we evaluated the successive 
integrations via the trapezoid rule until the solutions converged. The actual solution 
defining the overall process properties is given by the matrices at the end of the crystal 
Ua{z = f Ui{z = f Va{z = ^,uj,uj"), and Vtiz = f 

We checked the accuracy of the result in a variety of ways: At first we evaluated 
the canonical transformation conditions in (|B.1[) . (jB.2p . (jB.3p . and (|B.4[) . For example 
in the case of (jB.2p we evaluated: 

duj'Uaiz = -,UJ,Uj')Vb{z = -,Uj",Uj') 



= D('^ff\z^^,uj,uj') (F.l) 

and determined the distance of D^'^^^^^z — ^,uj,uj') from the expected zero matrix 
and consequently the error in the solution via: 

error = - - - 2 ' ' ' /--p 2^ 

0.5 [/ duj J duj'\Va{z = ^,uj,Lu')\ + / dw / duj'\Ubiz = ■|,a;,a;')|] ' 

In all presented cases the obtained error was below 0.00019. 

We also checked the numerical Schmidt decompositions of Ua{z,uj,uj"), 
Ub{z,uj,uj"), Va{z,uj,uj"), and Vb{z,uj,uj") to verify the mode properties de- 
rived in [Appendix B[ During this process the decompositions of Ua{z,u!,u!"), 
Ub{z,ui,uj") showed numerical issues, these however could be resolved by de- 
composing Ua{z,u},uj")Ul{z,u},u}"), Ul{z,u:,u:")Ua{z,u},uj"), Ub{z,uj,uj")Ul{z,uj,uj"), 
U^{z,u!,u!")Ub{z,u!,u!") iustcad. The obtained modes from these four matrices pro- 
vided a much improved stability especially in the high gain regime. Using these 
Schmidt modes we verified that the obtained Schmidt values of the U and V matrices 
behaved like cosh(rfe)^ — sinh(rfe)^ = 1 with errors below 0.00006. We also asserted that 
the decompositions yielded the functions ipkyipk, 4'k, £,k with symmetries as detailed in 
(|B.5p . which where fulfilled within numerical accuracy. 

The program code, published on our website, is able to directly create the 
investigated parametric down-conversion processes and also performs all mentioned 
tests. 
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